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Abstract: We propose an extension of the Fokker-Planck model of the Boltzmann equation to 
get a correct Prandtl number in the Compressible Navier-Stokes asymptotics. This is obtained 
by replacing the diffusion coefficient (which is the equilibrium temperature) by a non diagonal 
temperature tensor, like the Ellipsoidal-Statistical model (ES) is obtained from the Bathnagar- 
Gross-Krook model (BGK) of the Boltzmann equation. Our model is proved to satisfy the properties 
of conservation and a H-theorem. A Chapman-Enskog analysis and two numerical tests show that 
a correct Prandtl number of | can be obtained. 
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1 Introduction 

Numerical simulations of rarefied gas flows are fundamental tools to study the behavior of a gas in 
a system which length is of same order of magnitude as the mean free path of the gas molecules. 
For instance, these simulations are used in aerodynamics to estimate the heat flux at the wall 
of a re-entry space vehicle at high altitudes; numerical simulations are also used to estimate the 
attenuation of a micro-accelerometer by the surrounding gas in a micro-electro-mechanical system; 
a last example is when one wants to estimate the pumping speed or compression rate of a turbo- 
molecular pump. 

Depending on the problem, the flow can be in non-collisional, rarefied, or transitional regime, 
and several simulation methods exist. The most common method is the Direct Simulation Monte 
Carlo method (DSMC) proposed by Bird [3j: it is a stochastic method that simulate the flow 
with macro particles that mimic the real particles with transport and collisions. It works well for 
rarefied flows, but can be less efficient what the flow is close to continuous regime (even if there are 
other version of DSMC that are designed to work well in this regime: see for instance HBJS). For 
non-collisional regimes, the Test Particle Monte Carlo (TPMC) method is very efficient and is a 
common tool for vacuum pumps 0. For transition regimes, there exist mature deterministic solvers 
based on a direct approximation of the Boltzmann equation that can be used mmmmmm- 
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We refer to m or further references on these solvers. Among these solvers, some of them are based 
on model kinetic equations (BGK, ES-BGK, Shakhov): while these models keep the accuracy of the 
Boltzmann equation in the transition regime, they are much simpler to be discretized and give very 
efficient solvers PH El [HIM. These models simply replace the Boltzmann collision operator by a 
relaxation operator toward some equilibrium distribution, but still retain the elementary properties 
of the Boltzmann operator (collisional invariants, Maxwellian equilibrium). Two of them (BGK and 
ES-BGK) also satisfy the important entropy dissipation property (the so-called “H-theorem”) |2j. 

Very recently, in a series of paper mmm\m Jenny et al. proposed a very different and 
innovative approach: they proposed to use a different model equation known as the Fokker-Planck 
model to design a rarefied flow solver. In this model, the collisions are taken into account by a 
diffusion process in the velocity space. Like the model equations mentioned above, this equation 
also satisfy the main properties of the Boltzmann equation, even the H-theorem. Instead of using 
a direct discretization of this equation, the authors used the equivalent stochastic interpretation of 
this equation (the Langevin equations for the position and velocity of particles) that are discretized 
by a standard stochastic ordinary differential equation numerical scheme. This approach turned 
out to be very efficient, in particular in the transition regime, since it is shown to be insensitive to 
the number of simulated particles, as opposed to the standard DSMC. However, this Fokker-Planck 
model is parametrized by a single parameter, like the BGK model, and hence cannot give the 
correct transport coefficients in near equilibrium regimes: this is often stated by showing that the 
Prandtl number-which the ration of the viscosity to the heat transfer coefficient-has an incorrect 
value. The authors have proposed a modified model that allow to fit the correct value of the 
Prandtl number [12], and then have extended it to more complex flows (multi-species [10] and 
diatomic mi)- However, even if the results obtained with this approach seem very accurate, it is 
not clear that these models still satisfy the H-theorem. 

In this paper, we propose another kind of modification of the Fokker-Planck equation to get 
a correct Prandtl number: roughly speaking, the diffusion coefficient (which is the equilibrium 
temperature) is replaced by a non diagonal temperature tensor. This approach it is closely related 
to the way the BGK model is extended to the ES-BGK model, and we call our model the ES-Fokker 
Planck (ES-FP) model. Then, we are able to prove that this model satisfies the H-theorem. For 
illustration, we also show numerical experiments that confirm our analysis, for a space homogeneous 
problem. 

The outline of this paper is the following. In section [2j our model and is properties are proved, 
including the H-theorem. A Chapman-Enskog analysis is made in section [3] to compute the trans¬ 
port coefficient obtained with this model. In section [4[ we show that our model can include a 
variable parameter that allow to fit the correct Prandtl number. Our numerical tests are shown in 
section [5] and conclusions and perspectives are presented in section [6} 

2 The Fokker-Planck model 

2.1 The Boltzmann equation 

We consider a gas described by the mass density of particles f(t,x,v) that at time t have the 
position x and the velocity v (note that both position x and velocity v are scalar). The corre¬ 
sponding macroscopic quantities are ( p,pu,E ) = ((l,u, ^\v\ 2 )f), where p, pu , and E are the mass, 
momentum, and energy densities, and (4>) = J R3 (j>(v) dv for any velocity dependent function. The 
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temperature T of the gas is defined by relation E = \p\u\ 2 + | pRT , where R is the gas constant, 
and the pressure is p = pRT. 

The evolution of the gas is governed by the following Boltzmann equation 

dtf + v ■ V x f = Q(f, /), (1) 


with 


and 


Q(f, f)(v) = 


1«*£K 3 Ja£S 2 


- f(v*)f(v)) r 2 |u - 14 1 dcrdv*, 


, v + w* \v — vJ 
v = —x-1--- 


, V + u* 
V * = 


2 

\v — V* 


a. 
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It is well known that this operator conserves the mass, momentum, and energy, and that the 
local entropy H(f) = (/log/) is locally non-increasing. This means that the effect of this operator 
is to make the distribution / relax towards its own local Maxwellian distribution, which is defined 

by 


M(f)(v) = 


p 


■ exp 


\v — u[ 

2 RT 


(2nRT) 3/2 

For the sequel, it is useful to define another macroscopic quantity, which is not conserved: the 
temperature tensor, defined by 


0 := - ((v — u) <g> (y — u)f) . 
P 


(2) 


In an equilibrium state (that is to say when / = M(/)), © reduces to the isotropic tensor RTI. 


2.2 The ES-Fokker Planck model 

The standard Fokker-Planck model for the Boltzmann equation is 

d t f + v ■ VJ = U/ v • ((v - u)f + RTV v f ), (3) 

see [5j. Our model is obtained in the same spirit as the ES model is obtained from a modification 
of the BGK equation: the temperature that appears in ([3]), as a diffusion coefficient, is replaced by 
a tensor II so that we obtain 

dtf + v-V x f = D(f), (4) 

where the collision operator is defined by 

D(f) = -V v - ((v-u)f + UV v f), (5) 

T 

where r is a relaxation time, and II is a convex combination between the temperature tensor 0 
and its equilibrium value RTI, that is to say: 

n = (1 - is)RTI + u@, (6) 

with v a parameter. According to El,n is symmetric positive definite if v e] — 1]. In fact, this 

condition is too restrictive and we have the following result. 
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Proposition 2.1 (Condition of definite positiveness of II). 

The tensor II is symmetric positive definite for every tensor 0 if, and only if, 


RT 


Xmax ~ RT 


< V < 


RT 


RT - X r 


( 7 ) 


where A max and \ m in are the (positive) maximum and minimum eigenvalues ofO. Moreover H is 
positive definite independently of the eigenvalues of 0 as long as: 

~ \ < v < 1, (8) 

Proof. Since 0 is symmetric and positive definite, it is diagonalizable and its eigenvalues are pos¬ 
itive. Moreover, because of the definition of II, both 0 and II are diagonalizable in the same 
orthogonal basis. Then there exists an invertible orthogonal matrix P such that 


and 


PUP T = 


Xi 0 0 

POP 1 = I 0 a 2 0 
0 0 a 3 


(1 — u)RT + vX\ 0 

0 (1 — v)RT + VX 2 


0 0 (1 — v)RT + 1 /A 3 

Clearly II is positive definite if and only if (1 — v)RT + nXi > 0 for every i, which leads to: 

RT RT 

<u < 


A max ~ RT RT — Xmin 

which is ([7]). At worse A m * n is equal to zero so that 

RT 

v < - = 1. 

RT - 0 

At worse X max is equal to Trace(Q ) = 3 RT so that 

RT 1 

v > - 


3 RT - RT 


On the contrary when X max tends to RT near Maxwellian equilibrium (note that X m in also tends 
to T in this case), v becomes unbounded. The last two inequalities provide inequalities ([8]). 

□ 

Proposition 2.2. The operator D has two other equivalent formulations: 

D(f) = ^V v -(nG(f)V v J—Y (9) 

and 

DU) = \v v - (n/v„ log (^j-^ ) , (io) 
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where G(f) is the anisotropic Gaussian defined by 


G(f) = 


exp 


(v — u)n — u) 


( 11 ) 


Y / det(27rII) 
which has the same 5 first moments as f 

((M, l\v\ 2 )G(f)) = (p, pu, E), 
and has the temperature tensor ((v — u) (g> (v — u)G(f)} = II. 

Proof. Relations (J 9 ]) and (10) are simple consequences of the relation \7 v G(f ) = —II" 1 (v-u)G(f). 
The moments of G(f) are given by standard Gaussian integrals. □ 


Now, we state that D has the same conservation and entropy properties as the Boltzmann 
collision operator Q. 

Proposition 2.3. We assume v satisfies 0 and that v < 1. The operator D conserves the mass, 
momentum, and energy: 

((l,v,l\v\ 2 )D(f))=0, (12) 

it satisfies the dissipation of the entropy: 


(D(f) log /) < 0, 


and we have the equilibrium property: 


D(f) = 0& f = G(f) & f = M(f). 

Proof. For the conservation property, we take any function <j>[v) and we compute 

{Dim = -\ {{{v - u)f + nv„/) • V v fi). 

With <j>(v) = 1, we directly get ( D(f)) = 0. With <f>{v) = u,, we get 

(D(f)vi) = uf)f + II ijd Vj f) = 0, 

since (vif) = pui = p(vif) and {d Vj f) = 0. For fi{v) = \\v\ 2 , we get 

(D(f)vi) = {{{vi - Ui)f + II i,jd Vj f)vi) = -^(3 pRT - pTrace(II)) = 0, 

since Q shows that Trace(II) = (1 — u)3RT + i/Ti’ace(0) and ([ 2 ]) implies TVace(0) = 3 RT. 


For the entropy we need some inequalities to proceed. As noticed in the proof of proposition [2TTJ 
© and II are symmetric and can be diagonalized in the same orthogonal basis. Like it is done in 
jl|, we will now suppose that we work in this basis so that: 


/ Ai 0 0 \ 

0 = 0 A 2 0 , 

V 0 0 As / 
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and 

/ (l — v)RT + vX\ 0 0 

n = I o (i — is)rt + v \2 o 

\ 0 0 (1 — v)RT + 

For the following it is important to remind that ^ = 3-RT and that all eigenvalues A i are 

strictly positive. 

Since we want to prove that (D(f)lnf) < 0, we compute (D(f)ln f) with an integration by 
parts to get 


= - E (1 - ^ T + ^ ((fa -«,)/ + W) (% 


1=1 

3 


Eb 


2=1 
3 


A,; 


(1 — v)RT + u\ 
Xi 


((Vi ~ Ui) f 


dif 

f 


— v)RT + vXi / ! 

z_ -3^-^((w* -“*)/ + Xidif) y-j- + x 


dif (yi - Ui) (vi - Ui) 


E( 1 - (i-")f T + FA, l<(( ^-“. )/ 


i=l 

= A + B + C, 


dif 

f 


A < 


with 


^ (1 - v)RT + u\i / a t , x a t\ f di f , ( u < ~ “*) 

A = - 2_^ - — - ( {(Vi -Ui)f + Xidif) I — + -—- 

i= l \ \ J 1 


* = -E 


2=1 

3 


(1 — u)RT + vX i 

~Xi 


((vi - Ui) f + XAf) - 


(Vj ~ Ui) 
Xi 


2=1 


dif 


g = -E( 1 - (1 ^ ) > T+1 ' A, K(fa-"^lt 


(13) 

(14) 

(15) 


First, note that A is clearly negative since 


,, a f , \ o i dif , (v>i ~ Ui) 
((Vi -Ui)f + Xidif) ( — H-—- 


„, , dif (vi — uA 

= /Ai( -J-+ 


and ((l — v)RT+uXi) and A* are positive (due to the assumption on u and proposition |2.1[ ), and / is 
positive too. Then, B is equal to 0 since by integration by parts, we get (^{{vi — uf)f + Xidif) 

~Ti (i v i ~ u i) 2 f) + (/)> which is 0 since ((vi - u^ 2 f) = pX t . 
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The last term C satisfies 


- -E > 


2=1 
3 


(1 — v)RT + 

X 


1 1 i - 11 ) / ( —f 




2=1 


A,; 


E 1 

i=l v 
3 

Ed 


(1 — u)RT + z/Aj 

X 


(/> 


- 10 


2=1 


-f)<» 


•-r? »--)»■ 


2=1 


The Jensen inequality gives ^ X — 1 = Wr' which shows that C is non 

positive (since v < 1), and hence concludes the proof of the entropy inequality. 

In the general case, 0 and II are not diagonal. However, the same analysis can still be done: 
it is sufficient to use the change of variables v' = P T v, where P is the matrix of the orthonormal 
basis of eigenvectors in which 0 and n are diagonal. Indeed, let us define = f(v), then 

V v f(v) = FVJ'W), and 

[ (1, v 1 , {y' - u!) <g> (y' - u'))f'(v') dv' = (p, pii, p®'), 


where v! = P T u, and where &' = P T ®P is diagonal. We also define n 7 = P T HP which is 
also diagonal ad we have n 7 = (1 — v)RTI + v&. Then it is easy to find that we again have 


( D{f) In(/)) = A + B + C, but now with A, B, and C that have to be defined by (13) (15) with 
prime variables f ,v', and where (} denotes integration with respect to v'. Then we are back to the 
diagonal case and the previous proof is still valid. 

For the equilibrium property, we use ([£]) to obtain 


D(f) 


f 

G(f ) 


= ( nV < 


/ 

'G(f) 


V, 




since n is symmetric positive definite. Moreover, this integral is zero if and only if = 0, which 

is equivalent to / = aG(f). Since / and G(f) have the same first 5 moments, this relation is true 
if and only if a = 1, and hence / = G(f). Finally, this last relation implies {(v — u) <8> (v — u)f) = 
((v — u) (g) (v — u)G(f )), that is to say 0 = n. Then using ([b]) gives 0 = RTI which implies 
G(f) = M(f) and hence / = M(f). Conversely, if / = M(/), it is obvious that G(f) = M(f) and 
D(f) = 0. ' □ 


3 Chapman-Enskog analysis 

In this section, we prove the following proposition. 
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Proposition 3.1. The solution of the kinetic model (21) satisfies, up to 0{e 2 ), the Navier-Stokes 
equations 


d t p + V • pu = 0, 

dt.pu + V • {pu (g) u) + Vp = —V • a, 
dtE + V • {E + p)u = —V • q — V • ( au ), 

where the shear stress tensor and the heat flux are given by 

2 

a = —p(Vu + ( Vu) T -V • u ), and q = —kV • T , 

3 

with the following values of the viscosity and heat transfer coefficients 

P = 1/7 v and K = \ t P R - 

2(1 — v) 6 

Moreover, the corresponding Prandtl number is 

■3 

Pr = 


2(1 — v) 


and e is the Knudsen number defined below. 


In the first step of the proof, we write the conservation laws derived from (21) and (12) 


dtp + Vz • pu = 0, 

d t pu + V x ■ {pu <g> u) + V x • £(/) = 0, 
d t E + V, • {Eu + Ti{f)u + q{f)) = 0, 


(16) 


(17) 


(18) 


(19) 


where £(/) and q{f) denote the stress tensor and the heat flux, defined by 

E{f) = {{v - u) <g> {v - u)f) q{f) = {\{v - u)\v - u\ 2 f) . (20) 


The Chapman-Enskog procedure consists in looking for an approximation of £(/) and q{f) up to 
order 0{e 2 ) as functions of p,u,T and their gradients. 

To do so, we now write our model in a non-dimensional form. Assume we have some reference 
values of length xy, pressure p*, and temperature T*. With these reference values, we can derive 
reference values for all the other quantities: mass density p* = velocity u* = \JRT *, time 

f* = x*/v*, distribution function /* = p*/(i?T*) 3 / 2 . We also assume we have a reference value for 
the relaxation time t*. By using the non-dimensional variables w' = w/w* (where w stands for any 
variables of the problem), our model can be written 


dtf + v- V : 


J = ~D(f ), 


( 21 ) 


where e = is the Knudsen number. Note that since we always work with the non-dimensional 
variables from now on, these variables are not written with the ’ in (21). 

Note that an important consequence of the use of these non-dimensional variables is that RT 
has to be replaced by T in every expressions given before. Namely, now II is defined by 


n = (1 — v)TI + i/0, 


(22) 









the temperature is now defined by 
and the Maxwellian of / now is 

M(f) = 


E = \p\u\ 2 + | pT, 


(27 tT) 


3/2 


exp 


\v — u I 
IT 


Now, it is standard to look for the deviation of / from its own local equilibrium, that is to say 
to set / = M(f)(l-\-sg). However, this requires the linearization of the collision operator D around 
M(f), which is not very easy. At the contrary, it will be shown that it is much simpler to look for 
the deviation of / from the anisotropic Gaussian G(f) defined in ©• Since it can easily be seen 
that M(f ) and G(f) are close up to 0(e) terms, this expansion is sufficient to get the Navier-Stokes 
equations. 

First, we give some approximation properties. 


Proposition 3.2. We write f as f = G(f)( 1 + eg). Then we have: 

{0-,v, \\v\ 2 )G(f)g) = 0. 

Moreover, if we assume that the deviation g is an 0(1) with respect to e, then we have 

U = TI + O(e) 

and 

G(f) = M(f) + 0(e). 


(23) 

(24) 

(25) 


Proof. We have already noticed that ((1, v, ^\v\ 2 )G(f)) = ((1, v, 3 M 2 )/) (see proposition 2 . 2 k, 
which gives (23). The other relations (24) and ( |25[ ) are obtained by simple Taylor expansions 
of ( 22 ) and ( 11 ) with respect to e. □ 

Now, we show that £(/) and q(f ) can be approximated up to 0(e 2 ) by using the previous 
decomposition. 


Proposition 3.3. We have 


and 


W)=pi + 


1 - V 


T(M(f)g) + 0(e 2 ), 


q(f) = eq(M(f)g) + 0(e 2 ). 

Proof. By using / = G(f)( 1 + eg), we have 

£(/) = £(G(/)) + eE(G(f)g) = pU + eT(G(f)g). 


(26) 

(27) 

(28) 


But since by definition £(/) = p0(/), then (22) also implies pH = (1 — v)pTI + vT, (/). Us¬ 
ing this relation in (28) gives £(/) = (1 — v)pTI + z/£(/) + eT(G(f)g) which yields ( pt| ). The 
approximation (27) is immediately deduced from the decomposition / = G(f)( 1 + eg). 

□ 
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Consequently, the Navier-Stokes equations can be obtained provided that £(M( f)g) and q(M ( f)g) 
can be approximated, up to O(e). This can be done by looking for an approximation of the deviation 
g itself: it is obtained by using the decomposition / = G(f)( 1 + eg) into (21), which gives 

d t G(f) + V ■ V x G(f) + 0(e) = - £ D{G{f){ 1 + eg)). 

First, we state what the expansion of D(G(f)( 1 + eg)) is. 

Proposition 3.4. We have 

D(G(f)( 1 + eg)) = e—M(f)Lg + 0(e 2 ), 
r 


where L is the linear operator 


Lg = 


M{f) 


V„ • (TM(f)'V v g). 


Proof. Formulation ([9]) of D(f) gives 


D{f) = -v„ • ( nG(/)v 


G(f)(l + g ff) 

’ G(f) 

A direct computation gives V v G ^G(f)" 9 ^ = £ ^^9- Using this relation and I 
D(G(f)( 1 + £ g)) = -V„ • (TG(f)V v g) + 0(e 2 ). Using (25) gives the final resu 


(29) 


(30) 


24) into (30) gives 
t. □ 


This proposition shows us that g satisfies M(f)Lg = r(dtG(f) + v ■ G(f)) + O(e), and then 

T 


using (25) gives 


Lg = uu) [dtM{f) + v ' v * M(/)) + ° (e) - 


(31) 


We now have to solve this equation approximately to obtain an approximation of g. 

First, the right-hand side of (31) is expanded, and the time derivatives of p,u, and T are 
approximated up to O(e) by their space gradients by using (19). Then we get 

j^UtMU) + v • V X MU)) = A(V) • + B(V) : Vu + 0(e), 

where V = (v — u)/VT and 


A(V) = 


\V\ 


2 ~\)V and B(V) = V®V--\V\ 2 I. 


Then (31) reads 


Lg = r[ A(V) ■ ^ + B(V) : Vu ) + O(e). 


(32) 


This equation can be solved by showing that A and B are eigenvectors of L. Indeed, we have the 
following proposition. 
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Proposition 3.5. The components of A and B satisfy 


LAi = —3 Ai and LBj 1 = —2 Bi 


Proof. By using the change of variables V = (v — u)/VT, the dehnition (29) reduces to 

Lg = V v ■ (Mo(V)Vvg), (33) 

where Mq(V ) = —T-g- exp(— ^y-). The direct computation of LA{ and LBj j is easily obtained 
from (33) and is left to the reader. □ 


This property shows that we can look for an approximate solution of (32) as a linear combination 
of the components of A and B. We find 

9 = ~~- A(V) - T -Vu : B(V) + 0(e). 

Then we just have to insert this expression into H(M(f)g) and q(M(f)g) to get, after calculation 
of Gaussian integrals, 


S(M(/)s) = -y (Vu + (Vu) T - • ul) + 0(e) 

q (M(f)g) = -^VT + 0(e). 


Now we use proposition 3.3 and we come back to the dimensional variables to get 


£(/) = pi — g(V« + (Vu) T — -V • u ), and q = —kV • T, 

O 

with the following values of the viscosity and heat transfer coefficients 


I 1 = 


Tp 


2(1 -v)' 


and k = -rpR. 
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Using these relations into the conservation laws (19) proves that p,pu, and E solve the Navier- 
Stokes equations (16) up to 0(e 2 ) with transport coefficients given by (18) and that the Prandtl 
number is 

phR _ 3 

T “ 2(1 - is )' 


Pr = 


4 The ES-Fokker Planck model with a non constant v 

This short section establishes how we deal with the Prandtl number in all cases. 
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4.1 Some limits of the model with a constant v 


In the previous section, we have found that the Prandtl number obtained for the ES-Fokker Planck 
model is 


H 5R _ 3 

~k~Y ~ 2(1 - v)' 


(34) 


and can be adjusted to various values by choosing a corresponding value of the parameter v. 

Moreover, we have seen that the model is well defined (that is to say that the tensor II is 
positive definite for every /) if, and only if — ^ < v < 1. This last condition leads to the following 
limitations for the Prandtl number: 


1 < Pr < -poo, 


so that the correct Prandtl number for monoatomic gases (which is equal to |) cannot be obtained. 
In the next section, we show that this analysis, which is based on the inequality ([8]) is too restrictive, 
and that there is a simple way to adjust the correct Prandtl number. 


4.2 Recovering the good Prandtl number 


The previous analysis relies on inequality of proposition 2.1 that does not take into account 


the distribution / itself: this inequality ensures the positive definiteness of II independently of /. 
However, proposition 2.1 also indicates that n is positive definite for a given / if u satisfies Q. 
This inequality is less restrictive, since it depends on / via the temperature T and the extremal 
eigenvalues of 0. 

Now, our first idea is that u can be set to a non constant value (it may depend on time and 
space): it just has to lie in the interval [— A RT , 1], so that it satisfies the assumptions of 


proposition 2.3 The second idea is that the Prandtl number makes sense when the flow is close 
to the equilibrium, that is to say when / is close to its own local Maxwellian M(/): in such case, 
0 = RTI + 0(e), and all the eigenvalues of 0 are close to RT up to 0(e) terms, which implies 
A max = RT + Oe . Consequently, the value of v now lies in [— 1] which shows that v can take 

any arbitrary value between —oo and 1 when e is small enough, and in particular the value v = — | 
that gives the correct Prandtl number Pr = | can be used. 

In other words, by defining v as a non constant value that satisfies ([7]) and is lower than 1, 
we can adjust the correct Prandtl number provided that v can be set to — | near the equilibrium 
regime. 

This analysis suggests a very simple definition of v : we propose to use the smallest negative v 
such that: 

• n remains strictly definite positive, 

This leads to the following definition: 


v = max —, — 


RT 


4 A max ~ RT 


(35) 


Note that in most realistic cases u is equal to — |. Indeed, v ^ — | implies A max > 1-8 RT which 
can only happen in case of highly non equilibrium flow with strong directional non isotropy: such 
cases are very specific and are usually not observed in aerodynamical flows, for instance. 
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5 Numerical tests 


We have chosen to present a stochastic method to implement the model, since it is well suited 
to diffusion operators [15] (this was also used in [T3]). In this paper we only present results in 
homogeneous cases to understand how it works. Further studies will be done in inhomogeneous 
cases and comparisons with full Boltzmann equation and/or ESBGK model will be provided in 
forthcoming paper. In this section we solve 

d t f=-V v -{(v-u)f + UV v f). (36) 

r 

Our goal is to show numerically that we are able to capture the correct Prandtl number for 
monoatomic gases. This numerical illustration is based on the following remarks. First, the density, 
velocity, and temperature of / are constant in time: this is due to the conservation properties of 
the collision operator. Second, the heat flux q satisfies 


— q = - q 

at t 


so that 


q(t) = q( 0) exp- 


3 1 


Finally, the tensor 0 satisfies the relaxation equation 


= ( RTI-Q ), 


(37) 

(38) 

(39) 


which shows that it tends to RTI for large times, and hence v tends to the constant value —5/4 
(see (35)). Then for large times t > s, v can be assumed to be constant, and (39) can be solved to 
give 


0(t) = exp ( — 


2(1 — v){t — s) 


0(s) + ( 1 - exp ( - 


2(1 — u){t — s) 


RTI. 


Using equations ( [40] ) and (38), one gets 

ln(|0 M (t)-flr|) = 2(1 — v) 

MM 3 


1 

Pr 


(40) 


(41) 


for i = 1, 2, 3, so that we are able to recover the Prandtl number by looking at the long time values 
of 0(t) and q(t). 


5.1 A stochastic numerical method 

We use a DSMC method to solve the problem. For homogeneous cases, the probability density 
function is approximated with N numerical particles so that 

N 

f(t,v) ~ a^UiSv^t), 
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where Ui is the numerical weight of numerical particle i and Sypt) is the Dirac function at the 
particle velocity Vi(t). Moreover a is defined through the constant density p = f R3 f(t,v)dv by 

N 

p = a "jP U) t . 

i 

In the test cases presented here, all numerical weights Ui are equal, and we just have to de¬ 
fine the dynamics of the numerical particles. To solve diffusive problems, it is well-known that 
using Brownian motion is a good way to proceed (see m) : the corresponding stochastic ordinary 
differential equation is called the Ornstein-Uhlenbeck process that reads 

dVi(t) = -- ( Vi(t ) - u) + AdB{t) (42) 

T 

for each 1 < i < N. The quantity dB(t) is a three dimensional Brownian process. The matrix A 
has to satisfy AA 1 = II. Several choices are available for A. The obvious one would be to use the 
square root of II (which is a positive definite matrix): this requires to compute the eigenvectors of 
II and may lead to expensive computations. We find it simpler to use the Choleski decomposition 
because of the simplicity of the algorithm. For the time discretization we use a backward Euler 
method. The complete algorithm is the following: 

1. Approximate the initial data f(0,v) by w<5y.( t ), where N is the number of numerical 

particles and uj a constant numerical weight. The velocities are chosen randomly according 
to the particle density function /(0,.). 

2. Compute the tensor 0 

3. Compute the three real eigenvalues the positive definite tensor 0 with Cardan’s formula. 

4. Compute the Choleski factorization II = A T A of II. 

5. For i from 1 to N, advance the velocity V- 1 through the process: 



where At = t n+l —t n , Bi,B 2 ,B^ are random numbers chosen through a standard normal 
law. Since the scheme is explicit we enforce ^ < 0.1 to ensure stability: this leads to around 
0 p(, r time steps to reach the final time Tf. 

The scheme we propose here preserves mass but does not preserve momentum and energy, like 
most DSMC methods. However, these quantities are preserved in a statistical way. Moreover at 
the end of each time step, the distribution is renormalized to keep the mean velocity constant and 
to decrease the error in 0. 
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5.2 Numerical results 


We present two different test cases: one for which the correction of v (section 5.2.1) is not activated 
and one for which it is activated (section 5.2.2). For both cases, the characteristic time r is set to 
one, which is sufficient fore the illustration given here. We also set R = 1: in other words, we work 
in non-dimensional variables here. 


5.2.1 First test case 


We use 1 million particles. We choose three independent laws for the three components of velocity 
of the numerical particles: 

• the first component is equal to 100s 4 — 20 where s follows a uniform law between [0,1], 

• the second and third components of the velocity follow a uniform law between [—50, 50], 


The choice for the first component seems a little bit strange but we need to have a non zero heat 
flux at the beginning of the computations to be able to capture a characteristic time of variation 
for q. However, we have chosen distributions whose variances are of the same order so that the v 
which is used all along this computation is equal to — The final time is set to 1. 

As expected, we observe the convergence of the directional temperatures (the diagonal 
elements of 0) towards the temperature and the relaxation towards the Maxwellian (figure [T]). 
Moreover, all along the computation, the correction of v is not activated and hence the Prandtl 
number defined by (34) is always equal to | (figure [ 2 J) . Finally, the Prandtl number is computed 
by using linear regression on the logarithm curves of 0(f) and q{t ) (figure cn), see the discussion at 

^ 71 = 0.6428, which is 


the beginning of section |5| we get a numerical Prandtl number Pr n = 
close to |. 


4.5001 


5.2.2 Second test case 

We use 100 000 particles. As before we choose three independent laws for the three components of 
velocity of the numerical particles: 

• the first component is equal to 10 000x 4 — 2000 where x follows a uniform law between [0,1], 

• the second component of velocity follows a uniform law between [—50,50], 

• the third component of velocity follows a uniform law between [—50,50]. 

We have only changed the first component in order to make active the correction for v by having 
most of the thermal agitation in one direction. All the parameters are the same as before except 
the final time which now is 0.5 (we only want to see the change of regime for u). 

At time 0.5 we are not at equilibrium which explains why the distribution of the first component 
of velocity is not yet a Gaussian (figure[3]). We can still observe the exponential convergence towards 
the mean temperature of the diagonal components (figure [3]). When one looks more precisely at the 
behavior of v and the Prandtl number, one can see that the correction of v is activated so that it 
varies from —0.5 to —1.25 and consequently the Prandtl number varies from 1 to | (figure [ 4 ]) . One 
may also see that while log |g| is still linear , log |Tn — T| shows a slightly curved profile (figure [4]). 
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This is confirmed by the linear fitting of the logarithms (figure [5]). However, the numerical Prandtl 
number computed with the slope of this linearly fitted lines is then Pr n = = 0.717 which 

is not so far from This shows that activating the correction of the parameter v gives correct 
results. 

One may observe that getting Pr = 1 at the beginning of the computation is not correct. 
However, it has to be noticed that this happens because the flow is very far from equilibrium for 
small times: in such a regime, the Prandtl number has no clear physical interpretation, since the 
Chapman-Enskog is not valid, and hence the transport coefficients are not defined. 

6 Conclusion 

In this paper, we have proposed a modified Fokker-Planck model of the Boltzmann equation that we 
call ES-FP model. This model satisfies the properties of conservation of the Boltzmann equation, 
and has been defined so as to allow for correct transport coefficients. This has been illustrated by 
numerical tests for homogeneous problems. Moreover, it has been proved that the ES-FP model 
satisfies the H-theorem. 

The construction and analysis of ES-FP models for more complex gases (like polyatomic or 
multi-species), as well as inhomogeneous numerical simulations, will be presented in a forthcoming 
paper. Based on our proof of the H-theorem, we believe this theorem should still be true for such 
extended models. 

Acknowledgments. This study has been carried out in the frame of “the Investments for the 
future” Programme IdEx Bordeaux - CPU (ANR-10-IDEX-03-02). 
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Figure 1: Left: time evolution of the diagonal components T \\, T 22 , T 33 of the tensor 0 and of its 
trace T. Right: histogram of the first component of velocity at final time t = 1. 




Figure 2: Left: time evolution of v (defined by (35)) and Pr (defined by (34)). Right: time history 
of log |T — Till and log \q\. 
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Figure 3: Left: time evolution of the diagonal components T\\. T 22 , T 33 of the tensor 0 and of its 
trace T. Right: histogram of the first component of the velocity at time t = 0.5. 



Figure 4: Left: v and Pr along time. Right: convergence of the logarithms of |T — Tn| and q 
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Figure 5: Left: comparison between the curve of the logarithm of |T — T \\| ant its linear fitting. 
Right: comparison between the curve of the logarithm of |g| ant its linear fitting. 
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